Transcritical flow of a Bose-Einstein condensate through a penetrable barrier 

A.M. Leszczyszyn 1 ^ G.A. El 1 Yu.G. Gladush 2 ! and A.M. Kamchatnov^ 

1 Department of Mathematical Sciences, 
Loughborough University, Loughborough LE11 3TU, UK 
2 Institute of Spectroscopy, 
Russian Academy of Sciences, Troitsk, 
Moscow Region, 142190, Russia 



■ The problem of the transcritical flow of a Bose-Einstein condensate through a wide repulsive 
| penetrable barrier is studied analytically using the combination of the localized "hydraulic" solu- 
tion of the ID Gross-Pitaevskii equation and the solutions of the Whitham modulation equations 

■ describing the resolution of the upstream and downstream discontinuities through dispersive shocks. 
It is shown that within the physically reasonable range of parameters the downstream dispersive 
shock is attached to the barrier and effectively represents the train of very slow dark solitons, which 
can be observed in experiments. The rate of the soliton emission, the amplitudes of the solitons in 
the train and the drag force are determined in terms of the BEC oncoming flow velocity and the 
strength of the potential barrier. A good agreement with direct numerical solutions is demonstrated. 
Connection with recent experiments is discussed. 
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I. INTRODUCTION 

The problem of the fluid flow past an obstacle is of great importance for both normal fluids and superfluids. In 
dynamics of viscous compressible fluids it is closely related to the theory of viscous shocks (see, e.g., [l|, 0]). In the 
theory of surface water waves it has lead to the development of a detailed description of ship waves and corresponding 
\ drag forces (see, e.g., @). In superfluid dynamics, understanding of the nature of the critical velocity in the flow past 
l— ~~ 1 ■ an obstacle was crucially important for the development of the superfluidity theory [1, d, @ . Therefore it is natural 
| that this problem has attracted much interest in the context of dilute gaseous Bose-Einstein condensates (BECs) 
dynamics where formation of dispersive shocks and "ship waves" has been studied both theoretically @, H, OH HH, 
f-^. ' UH, [HI, 0, EH and experimentally [ill, E2i El ■ So far the main attention was paid to physically multi-dimensional 
\& , problems which in some cases could be asymptotically reduced to effectively one-dimensional (ID) models for the 
purposes of mathematical convenience. However, there exist essentially ID situations when the phenomena observed 
are specific to the ID case only (see, for instance, the discussion in [I3)- One of such prototypical problems arises in 
<r— ( ■ the description of the water flow over a bottom ridge [H, [h| [2(1 . In the BEC context, a similar problem arises in the 
study of the BEC flow induced by a penetrable barrier moving through an elongated BEC, as it has been recently 
studied in the experiment [2lj . The most characteristic phenomenon observed in this experiment is that a chain of 
dark solitons is generated for a finite interval, v_ < v < ill, of the barrier velocities v. The existence of certain 
^ , threshold velocity V- agrees with the main concepts of the superfluidity theory introduced by Landau [4] and with 
earlier experiments [22] on the appearance of a critical velocity for an "impenetrable" obstacle moving through the 
condensate at different velocities. Taking into account the fact that formation of vortices is impossible in a quasi-one- 
dimensional setting, one can easily deduce from the dispersion relation for linear Bogoliubov excitations in weakly 
non-ideal Bose gas that the critical velocity u_ coincides with the sound velocity, and this conclusion was confirmed 
by the experiment [22l |. However, nonlinear effects modify considerably the Landau criterion based on the idea of the 
generation of linear excitations at velocities greater than the critical one. For example, in the case of wide smooth 
potentials used in [5l[ one has to take into account nonlinear effects which make possible the generation of nonlinear 
excitations (dark solitons) and which reduce the critical velocity V- to a subsonic value [23| . It is remarkable that 
this nonlinear mechanism is effective for a finite interval of the barrier velocities only. As was first noticed in [24| . 
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FIG. 1: Sketch of experimental setup in Ref. [2jJ. The potential barrier created by the laser beam is moved through the BEC 
along the —x direction. The trap tightly confines the BEC in the radial y, z directions. 



there exist such special forms of the barrier potential that no radiation is generated at supersonic velocities greater 
than the second critical velocity v+, and this result was confirmed numerically for the potentials of a more general 
form. As was noted in [25], while generically some linear radiation still exists for v > v+, its amplitude decreases 
exponentially with the growth of the ratio of the obstacle size to the healing length so that "broad and smooth" 
potentials can always be considered as "radiationless" outside the interval [v_,t>+]. These results were also confirmed 
by the numerical simulations of waves generated by a supersonic motion of a repulsive rectangular potential (26j 
and oscillatory supersonic motion of the Gaussian potential [2?| . Thus, for smooth wide penetrable potentials the 
production of solitons is effective only in a finite range of the barrier velocities. In the present paper, motivated by 
the results of the experiment [2l[ and the theoretical setting of [HI, [2(| [23| , we develop full asymptotic theory of the 
transcritical BEC flow for the case of broad repulsive potentials. 

Although the experiment [2l[ was performed with a dense BEC for which the radial motion of the gas is essential, 
we shall confine ourselves here to the case of a rarefied gas with a frozen radial motion when the full Gross-Pitaevskii 
(GP) equation can be reduced to the ID nonlinear Schrodinger (NLS) equation (see, e.g., [28|]). This limiting case 
reproduces all the main characteristic features of the phenomenon, it admits exhaustive analytical treatment, and 
such a setup can be quite suggestive for future experiments. 

II. FLOW OF A BOSE- EINSTEIN CONDENSATE THROUGH A PENETRABLE BARRIER 

A. Mathematical model 

Engels and Atherton [2~H considered a wide penetrable barrier moving with constant speed through an elongated 
BEC confined to a cigar-shaped trap, as it is illustrated schematically in Fig. 1. In the experiment, the role of the 
barrier was played by a repulsive potential created by the laser beam which is swept through the BEC with the 
velocity v along the —x direction. If the condensate is rarefied enough, then its dynamics is described by the forced 
NLS equation which in standard non-dimensional units has the form 



iipt + hipxx - = V{x + vt)il>, 



(1) 



where V(x + vt) > is the moving potential barrier. It is convenient to transform this equation by means of the 
substitution 



ip(x, t) = \J p(x, t) exp ( i / u(x',t)dx / 



(2) 



to a hydrodynamic-like form and pass to the reference frame moving with the velocity —v by introducing x' 
x + vt, u' = U + V. 



li' t + U'u' x , + p x , + [ 



4P ) x , 



pt + {pu') x > = 0, 

+ 14,(0 = 0. 



(3) 
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Here p and u are the condensate density and velocity, respectively. It is supposed that in this reference frame and in 
our non-dimensional units the flow at infinity satisfies the conditions 

p — ► 1, u' —>■ v as \x\ — > oo, (4) 

that is the length of the elongated condensate is assumed to be much longer than the size of the wave structures 
generated in the flow. Thus, equations ©-(HI) represent an idealized mathematical model for the description of BEC 
dynamics in the configuration of our interest. 

The study of the problem ©-Q was initiated in Ref. [23| where the time-independent flows induced by a slowly 
varying in space barrier moving with subcritical velocity were analyzed. Unsteady supercritical flows were considered 
in [23| numerically and it was noticed there that the results are reminiscent of those for the flow of a stratified fluid 
over a broad localized topography [H, H3] • In @, Hil the shallow- water flow past topography problem was studied 
in the framework of the forced Korteweg-de Vries (KdV) equation, and this approach provides a clue to solving our 
NLS equation problem j3])-([4]). The first step in this direction is the study of the stationary solutions of Eqs. ©-(H]) 
in the case of a wide and smooth potential V(x) so that the dispersion terms in ([3]) can be neglected and one can 
take advantage of the so-called "hydraulic" or "hydrostatic" approximation. 



B. Hydraulic solution 



1. General construction 



Let the potential V(x) be localized in the interval (—1,1), where I 1 (i.e. the spatial range of the potential in 
dimensional units is supposed to be much greater that the healing length of the BEC). Then one can make use of 
the hydraulic approximation by assuming that the characteristic length at which all variables change has the order 
of magnitude of /. Thus one can neglect the terms with the higher derivatives in (|3|) and consider the stationary 
solutions described by the system 

(pu) x = 0, 
uu x + p x + V x {x) = 0, 

where we have omitted primes for convenience of the notation. Actually, these are equations for the stationary 
hydraulic flow in the reference frame with the barrier at rest (we note that the flow direction in our paper is opposite 
to that in 23]). Both equations ([5]) are readily integrated once to give 

pu = v, \u 2 + p + V{x) = \v 2 + 1 , (6) 

where the integration constants are found from the boundary conditions J4]). Equations ([6]) thus define a stationary 
flow that connects smoothly to p = 1, u = v at both infinities. Eliminating p we obtain an implicit equation for the 
flow velocity u as a function of x, 

V(x)=F(u), where F(u) = \(v 2 - u 2 ) - - + 1. (7) 

u 

For ([7]) to have a solution defined for all x one should require that 

V m < \v 2 - \V 2 ' Z + 1 , (8) 

where V m = maxV(i). Indeed, the function F(u) varies within the interval — oo < F(u) < p(v), where p(v) — 
max{F(u)} = jv 2 — |i> 2 / 3 + 1. So for the equation |7]) to have a real solution for all x the interval [0, Vm\ must lie 
within the range of the function F(u) (the condition © in the present superfluid context was obtained in [23, [25| but 
can also be found in earlier studies on shallow- water flows past topography — see, for instance [lj|). 

Inequality ([8]) defines the subcritical v < u_ and supercritical v > v + regimes where v± are the roots of the equation 
(see Fig. 2) 

p{v) - \v 2 - ft- 2 / 3 + 1 = V m (9) 

and we have V- = for V m > 1. If V m <^ 1, then v± must be close to unity so one can easily find that up to the 

1/2 

second order in the small parameter V m the critical velocities are equal to 



v ± *l±\l— + — . (10) 
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FIG. 2: Plot of the function p(v); for V m = 0.5 the critical values are equal to u_ « 0.2, v + ^ 1.9. 
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FIG. 3: Plots of the fluid velocity u(x) in the hydraulic solution (|lip for the potential barrier (|13l) for which v~ « 0.2, v+ ~ 1.9; 
(a) corresponds to subcritical velocity v = 0.1 and (b) to supercritical velocity v = 2.0. The physical branches are shown by 
solid lines and nonphysical by dashed lines. 



In the transcritical regime v_ < v < v + condition ([8]) does not hold so it is natural first to look closer at what happens 
with the hydraulic solution when v approaches its boundaries from outside. We now look at the dependence of the 
flow velocity u{x) and the density p{x) on the space coordinate x determined by equations ©, (O which can be 
re-written as 



\u 2 + ^ + V(x) = \v 2 + 1 (11) 



and 



2p 2 



+ p+V{x) = ±i/ + l, (12) 



respectively. Formally, the solution to (fTTj) (or (fT2"|0 has two branches only one of which satisfies the necessary 
boundary conditions ([4} so just this branch should be considered as the physical one. These two branches of the 
solution for the flow velocity u(x) are plotted in Fig. 3 for the case of the potential 

v & = It i \ (i3) 

cosh(x/(T) 

with V m — 0.5, a — 2, for which Eq. ([9]) yields w_ « 0.2, v+ « 1.9. We see that the subcritical (v < V-) and 
supercritical (v > v + ) hydraulic solutions have similar structure but are characterized by the "exchanged" relative 
positions of the physical and nonphysical branches. 

In Fig. 4 we have plotted two series of hydraulic solutions for different subcritical and supercritical values of the 
velocity v. As one can see, the physical and nonphysical families are separated by two separatrices. When one 
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FIG. 4: a) subcritical hydraulic solutions for u at different values of v < v- (thin solid lines) and a separatrix solution at v = v- 
(thick solid line); b) supercritical hydraulic solution at different values of v > «_ (thin solid lines) and separatrix solution at 
v — (thick solid line). Short-dashed lines: non-physical branches. Long-dashed line: non-physical separatrix. 



approaches v = i?_ from below (Fig. 4a), the physical branch becomes more and more pointed up at the center of 
the barrier (x = 0), and when v = it bifurcates into the separatrix line. This separatrix solution satisfies the 
necessary boundary condition v = V— asu — oo (effectively at x — —I) but it fails to satisfy the same equilibrium 
condition at x — > +oo so one can expect that this bifurcation will be accompanied by the generation of an unsteady 
flow downstream. 

Indeed, the analysis of the near-critical NLS flow through the delta-functional potential V(x) in (23|, [3(| shows 
that for some v = v cr (V m ) the flow loses its stability through the saddle-node bifurcation resulting in the generation 
of dark solitons for v > v cr . The numerical simulations in [23j | also suggest that one can expect a qualitatively similar 
scenario with the soliton train generation for slowly varying potentials moving with the supercritical velocity. 

To describe this generation of solitons in the NLS flow past broad potential barrier quantitatively, we shall take 
advantage of the analytical construction proposed in 20], where the transcritical shallow- water flow past extended 
localized topography was considered in the framework of the forced KdV equation. The key in this construction is 
the assumption (confirmed numerically a posteriori) of the existence, for certain interval of the flow velocities, of the 
local steady transcritical transition over the forcing region, which is described by the relevant hydraulic solution. This 
solution does not satisfy the equilibrium boundary conditions outside the spatial range [— I, I] of the forcing, so one 
introduces discontinuities at x — ±Z, which are resolved into the equilibrium state with the aid of unsteady nonlinear 
wave trains (undular bores). By applying a similar assumption here for v = i?_ we consider the described above 
separatrix hydraulic solution as a local one defined on the interval [— I, I]. Then we get a discontinuity at x — I, which 
should be further resolved into the undisturbed state w = w,p = 1 as i -> oo via the dispersive shock wave. Thus, 
when the BEC flow velocity is equal to v = w_ , the formation of the dispersive shock downstream the potential barrier 
is expected. 

In a similar way, when we get closer to v = v + from above (Fig. 4b), the physical solution becomes pointed down, 
and at v = v + it bifurcates into the separatrix line which satisfies the boundary condition u — v + as x — > +oo 
(effectively at x = I) and a discontinuity u = u_ ^ 14. occurs at the left edge x = —I. Hence, in this case v = v + we 
expect the formation of the upstream dispersive shock. One should note that in the NLS dispersive hydrodynamics a 
general discontinuity resolves into a certain combination of two waves: dispersive shock(s) and/or rarefaction wave(s). 
The actual combination depends on the relation between the specific initial jump values for the density and velocity 
[3lL |32| . The closure conditions enabling one to single out the unique combination in our problem will be formulated 
in the next section. 

Thus, one can suggest that in the two limiting cases v = v± the full solution can be built of two parts: the steady 
hydraulic transition solution defined within the spatial range of the potential barrier and the unsteady dispersive 
shock/rarefaction wave combination at one of the sides of the barrier. Since the dispersive shock is generated essentially 
outside the spatial range of the potential, one can use the potential-free NLS equation for its description and take 
advantage of the theory developed in I3ll. 13211 on the basis of the original Gurevich-Pitaevskii approach using the 
Whitham method of slow modulations [33, [3J] . This modulation theory of the NLS dispersive shocks was recently 
applied to the description of dispersive shocks in freely expanding BECs in [1, [TTJ] . 

Generally we are interested in the flow corresponding to the transcritical region 



w_ < v < v+ 



(14) 
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so one can expect formation of both upstream and downstream discontinuities with their further dispersive resolution 
into the undisturbed states at x — > ±00. To be precise, we expect a subcritical jump upstream, a supercritical 
jump downstream and the exact criticality at the top of the potential. One should emphasize that the suggested 
description of the transcritical BEC flow as a combination of local hydraulic solution connected with the equilibrium 
state at infinities by unsteady dispersive shock(s) is an assumption which should be validated by the comparison of the 
analytical solution with direct numerical simulations for a range of the oncoming flow/potential barrier parameters. 
As we shall see, this assumption works very well when the oncoming flow velocity is not too close to the transcritical 
region boundaries v±(V m ). The reason for that is clear: if the speed v is close to one of its critical values v±(V m ), the 
characteristic relaxation time is expected to be very large (see [30] for the bifurcation scaling analysis in the case of 
the delta-functional potential) so the local transcritical steady flow does not establish in finite time. Nevertheless, it 
is instructive to perform the analysis, based on the above assumption, for the whole range of velocities in the interval 
[i>_,i>+] and then see how well this approximation works for different parameter values. 
We denote the values of the density and velocity at the boundaries x = ±1 of the hydraulic solution as 

p(—l) = p u , u(—l) = u 11 (upstream the barrier) (15) 

and 

p(l) = p d , u(l) = u d (downstream the barrier) . (16) 

We assume that the potential maximum is located at x = 0, i.e. V x (0) = 0. The transcritical hydraulic solution of 
our interest is then distinguished by two sets of conditions: 

(i) u x ^ at x = 0; (17) 



(ii) u u < u m , p u > p m and u d > u m , p d < p m , (18) 

where u m , p m are the values of u and p at x — 0. As we shall see the condition (i) actually coincides with the 
requirement that u m is equal to the local sound velocity at the top of the potential (exact criticality) . 
Now we take the solution of the hydraulic equations |5]) in the general form 

pu = ci, i« 2 + p + V(x) — c 2 , (19) 

01,2 being arbitrary constants, so that instead of ([7|) we obtain 

C2-^--=V(x). (20) 
2 u 

Differentiating (|20|) we obtain 

(u ~^)u x = V x . (21) 

Applying condition (fT7| we get 

d = t& . (22) 
Combining this relation with the first integral (fT9|) applied to the point x — yields 



.2 



which simply means that the local sound velocity at x — is ^/p m = u m . 
Substituting (|22|) , (|23|) into the second integral (fT9| gives the value of c 2 : 



(23) 



C2 = \ul l + V rn . (24) 

On the other hand, the integrals p^|) applied to the upstream and downstream boundaries x = ±Z where V{x) — > 
yield the relations 

p^ d = ul, l( u ^ d f+p^ d = lu 2 m + V m . (25) 
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Eliminating u m we get three equations 

p u u u = p d u d , l(u u ) 2 + p u = ±(u d ) 2 + p d , \{u u Y+p v --l{p u u u f/ s = V m (26) 

for four quantities Obviously, the last equation can be replaced by 

l(u d ) 2 + p d ~l(p d u d ) 2 / 3 = V mi (27) 

that is u u , p u and u d , p d satisfy the same set of equations. One more equation is needed for closing the system, and 
this will be obtained in the next subsection. 

Once u u,d and p u,d are found, the integration constants c\ and C2 in (|19p are expressed as 

d = p u4 u u ' d , c 2 = \{u u ' d ) 2 + p u ' d . (28) 

Hence, the transcritical hydraulic solution (fT9|) is given by 

+ + V{x) = U u u, d) 2 + p u,d^ (2Q) 

z U z 

( n u,d u.d \ 2 
9 - — j +p + V(x) = ±{u u - d ) 2 +p uA . (30) 

At the boundaries v = v± of the transcritical region, where either p u = 1, u u = V- or p = 1, u d = u + , equations 
(HHJ), ((201), reduce to ([XT]). ([T2]l respectively as it should be. 

i?. Closure conditions 

To get the closure condition for (|2"6")l we shall take advantage of the transition condition across the dispersive shock 
which is most conveniently formulated in terms of the Riemann invariants of the dispersionless hydrodynamic system 
associated with the NLS equation. This system is obtained from ([3]) by removing the dispersion and potential terms, 

pt + (pu) x = 0, 
u t + uu x + p x = , 

and is equivalent to the classical shallow-water equations. Introducing the Riemann invariants 

A± = i u ±Vp, (32) 

we represent (|31j) in the diagonal form 

^ + i(3A + + A_)^^0, ^ r + i(A ++ 3A_)_ = 0. (33) 

In the problem of the decay of an initial discontinuity in the NLS hydrodynamics the transition across the dispersive 
shock is characterized by the zero jump of one of the hydrodynamic Riemann invariants (|32p [3ll |32|. It should be 
stressed that this condition is not a small-amplitude approximation of the classical shock jump condition but rather is a 
non-perturbative consequence of the data transfer along characteristics of the associated Whitham modulation system 
describing dispersive shock region [3^ |. It will transpire later that in the present problem of the right-propagating 
BEC through a potential barrier at rest, the "conserving" invariant is A+ (this corresponds to the waves propagating 
to the left which might seem somewhat counter-intuitive with regard to the downstream dispersive shock) . 

The upstream dispersive shock is adjacent to the oncoming flow p = 1, u = v so the Riemann invariant condition 
yields the relationship 

\u u + - \v + 1 , (34) 



which closes the system (1261) . We also recall that we are interested in the solution satisfying condition (| 18[) which 
allows one to identify the roots of |26|) . (|34|) as upstream and downstream ones. Now all four values p u,d , u u,d are fixed 
in terms of V m , v, which, in particular, implies that the downstream discontinuity generally cannot be resolved by 
a single dispersive shock but requires an additional right-propagating rarefaction wave to adjust to the undisturbed 
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FIG. 5: Riemann invariants A± as functions of x in the hydraulic solution (|29[) . (|30|l for the potential (|13|) with V m ~ 0.5. Note: 



flow at +oo. While the introduction of this rarefaction wave does not present one with a serious technical problem, 
we shall see that in practical terms this additional complication is not actually necessary. 

First we note that, since the matching of the (external) hydrodynamic solution and the (internal) modulation 
solution describing dispersive shock is most naturally formulated in terms of Riemann invariants (see (3lL [32j for 
details) it is instructive to plot the transcritical hydraulic solution (|19p , (|22p , (|24p in terms of A± (the physical branch 
of the solution is selected using inequalities (|18p). The plot in Fig. 5 suggests that for weak enough potentials one can 
neglect the change of A+ across the barrier, which would remove the necessity of introducing an additional rarefaction 
wave. To justify this supposition we first eliminate p u ' from the second equation (|25p to represent it in the form 

^ + £-§•4-"- < 36 > 

Introducing the normalized quantities v± = u < u /u m we represent (1331) as 

9 2 

v%.-\ 3 = a, (36) 

where a — 2V m /uf n . Now let us suppose that V m 1, which implies u± ~ 1, u « u m » 1, cv<C 1. Then from ((36]) 
we get the expansion 

» ± = l±(^) 1/2 + f + c ±a ^ + ..., (37) 

-|/2 

i.e. the controlling small parameter is again e ~ V m (note that the coefficients c± in (|3T[) will not contribute to the 
result so we do not present them explicitly) . We now consider the boundary values of the Riemann invariant A+ (|32p 

i 3/2 



^ = ^ M + V^=^ + 7 -5^. 08) 



Again, normalizing, A u ' d = A"' d /u m , and expanding for small V m we obtain 



Now, taking into account that u m — 1 to leading order, we get 



(40) 
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Thus, for weak potentials the jump of the Riemami invariant A+ across the potential has the third order in the 
controlling small parameter. This result has certain analogy with the classical result from the shock wave theory 
which reads that the relevant Riemann invariant has just the third order jump across the weak shock [l], Q]. The 
coefficient (f/6) 3 / 2 before Vm" 2 suggests that one can neglect the jump of A+ even for the potentials of moderate 

1 /2 

strength (of course, provided that the coefficients for the successive powers of the small parameter V m in the 
expansion (|40[) are reasonably small). Indeed, for our potential (|13p with V m — 0.5 we have A+ d w 1.5 while according 
to (|40p (5 w 0.025 i.e. just about 1.7%, which is confirmed by the numerical transcritical hydraulic solution shown 
in Fig. 5. So, for weak to moderate potentials one can safely assume that the value of the Riemann invariant A + is 
preserved across the potential, which allows one to use an additional closure condition 

\u U + V7 F = \u d + yfjfi = \v + l. (41) 

Relation (|4ip is asymptotically consistent with exact conditions ([26]) . (|34|) and will be especially useful in our further 
consideration of the dispersive resolution of the downstream discontinuity. We also note that an immediate implication 
of (|4"l"j) is that both upstream and downstream dispersive shocks are "based" on the same family of characteristics 
corresponding to left-propagating hydrodynamic simple waves, which will be essential for the modulation solution in 
the subsequent sections. 



Weak potentials: explicit formulae 



Using asymptotic closure conditions (|41[) one can obtain simple approximate explicit expressions for p u,d , u u ' d in 
terms of v,V m , which will be useful later. We use (|4Tj) to eliminate p u - d from (|26|) to obtain a single equation for 

..u,d 



W = U 



w 

T 



w 



1 2/3 



Vrr, 



(42) 



This equation has two roots, the larger one corresponds to u d and the smaller one to u u (see (|18p). 

Equation (|4"2"|) can be solved approximately for V m <C 1. It is easy to see that for V m = this equation is satisfied 
if w = 1 + (v — w)/2 i.e. w = (v + 2)/3. In the next approximation we obtain 



l+g(«-l) 



2V„. 



2Vj77 



where we have used that v takes its values in the transcritical region (see (fT0|l ) 

1 — \ — — < v < 1 - 



Respectively, with the same accuracy we get from (|4T 



1 + -(w - 1) + 



2V„, 



2V m 



(43) 



(44) 



(45) 



Using (|43l) . (I45p we calculate the upstream and downstream values of the Riemann invariant A_ = — ^fp, which 
undergoes discontinuities at both edges of the transcritical hydraulic solution (see Fig. 5), 
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(v-1) 



2V m 



MO 



12 



(v-l) + 



2V„ 



(46) 



Obviously, A" < A^°, X_ > X°° , where Af° = v/2 — 1 is the value of A_ for the undisturbed flow at infinity. 



C. Resolution of downstream and upstream discontinuities 



Summarizing the results of the previous section, we have the following values for the Riemann invariants A± at the 
boundaries x = ±1 of the hydraulic transition and at x = ±oo: 



M-oo) = A + = A + = A +(+oo) = |« + 1 , 



(47) 
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A_(-oo) = A_(+oo) = \v - 1 , (48) 



Al = i u d -vV, X u _ = \u u -^ (49) 

From now on we consider the values A"' d (v, V m ) as known. Thus, there are upstream and downstream discontinuities 
in A_ and no discontinuities in A+. 

As was mentioned above, the upstream and downstream discontinuities are resolved through the generation of 
dispersive shock waves which can be described using modulated periodic solutions of the defocusing NLS equation 
(P) without the barrier potential V(x + vt). The potential term can be neglected because the shock resolution occurs 
essentially outside the potential range —l<x<l. The potential-free NLS equation is Galilean invariant and hence it 
preserves its form (up to inessential in our case phase factor in ^-function) after the transformation to the reference 
frame with the barrier at rest which is used in our calculations. 

The periodic solution of the NLS equation can be written in terms of the condensate density as (see, e.g., [3l|, HBj) 



p{x,t) = ±(A 4 - A 3 - A 2 + Ai) 2 + (A 4 - A 3 )(A 2 - Ai)sn 2 ( ^(\ A - A 2 )(A 3 - Ai) £, m). (50) 

Here 

( -'- m ' u =#- (51) 

U being the phase velocity and m the elliptic modulus, < m < 1. The condensate velocity u is expressed in terms 
of density as 

u(x,t) = U--^-, (52) 
p{x,t) 

where 

C = i(-A! - A 2 + A 3 + A 4 )(-A : + A 2 - A 3 + A 4 )(Ai - A 2 - A 3 + A 4 ) . (53) 
The wavelength is expressed in terms of A^'s as 

L= - ^ , (54) 

v/(A 4 - A 2 )(A 3 -Ai) 

K(m) being the complete elliptic integral of the first kind. 

The parameters Ai < A 2 < A 2 < A 4 change slowly through the dispersive shock and their evolution is governed by 
the Whitham modulation equations 

^ + ^(A)|^ = 0, i = l,2,3,4, (55) 

so Aj's are the Riemann invariants of the Whitham equations. Here Uj(Ai, . . . , A 4 ), i = 1,2, 3, 4, are the characteristic 
velocities which can be expressed in terms of the complete elliptic integrals of the first and the second kind [13, HH ■ 
One can use the following convenient formula for the computation of w,'s [35l |36| 

The waveform (1501) within the dispersive shock wave region gradually changes from the vanishing amplitude har- 
monic wave at one of the edges to a dark soliton at the opposite edge so that generally the modulus m runs over 
the whole range from to 1 . The Riemann invariants Xj of the Whitham equations are matched with the "external" 
hydrodynamic invariants X± at some free boundaries x^(t) defined by the conditions m — (harmonic edge) or m = 1 
(soliton edge) . The specific matching conditions and relative position of the harmonic and soliton edges depend on 
whether one considers the left- or right-propagating dispersive shock (see [HJ ). For the left-propagating case of our 
interest the matching conditions are 

A 2 = Ai , A 4 = A+, A 3 = A_ at x = x~(t); . . 

A 2 = A 3 , A 4 = A + , Xi — A_ at x — x + (t). ^ ' 
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FIG. 6: Qualitative behavior of the Riemann invariants Xj as functions of x in the full hydraulic/modulation solution. 



Since dispersive shocks expand with time, their widths at sufficiently large t are much greater than the range / of 
the barrier potential, so for t I one can assume the hydraulic solution to be asymptotically localized at x — and 
use the similarity solutions of the Whitham equations (|55ll with only one variable changing through each shock (3lj | 
(note that the formal condition of applicability of the Whitham method itself to the decay of a step problem is t 3> 1 

— see [H, Hl|)- From the matching conditions |57|) it is clear that the varying invariant should be A2. Thus, in each 
dispersive shock the invariants Ai, A3, A4 are constant and the invariant A2 varies according to the similarity solution 



x 

w 2 (Ai, A 2 ,A 3 ,A 4 ) = -, (58) 



where 



r\ \ \ \ \ iTi (A 3 - A 2 )(A 2 - Ai)if(m) 

u 2 Ai, A 2 ,A 3 ,A 4 ) = 5 >A + — rrr^-T — 7T \~\~Fpl T' 59 ) 

^ (A 3 - \i)K (m) - (A 3 - Ai)-E(m) 

E(m) being the complete elliptic integral of the second kind. The values of Ai, A3, A4 in (|59[) are fixed and the Riemann 
invariant A4 must have the same constant value 

A 4 = 5^ + 1 (60) 

in both downstream and upstream dispersive shocks since the value of A+ is transferred through the hydraulic 
transition (see (|47)l ). At the same time, the values of Ai and A3 are different upstream and downstream the barrier 
so we shall consider the downstream and upstream dispersive shocks separately. Before we proceed with the detailed 
analysis of the dispersive shocks it is instructive to get a picture of qualitative behavior of the Riemann invariants 
in the entire flow. For that, we use the matching conditions ([57]) and the inequality 8X2/ dx > following from the 
similarity modulation solution (f58|) and the general property dvi/dXi > of the characteristic velocities (l56l) . As a 
result, one arrives at the scheme of the behavior of the Riemann invariants sketched in Fig. 6. 

1. Downstream dispersive shock/ soliton train 

In the downstream dispersive shock we have (see the matching conditions (|57p and Fig. 6) 

A 1 =A~ = |-1, A3 = Al, A 4 = A^ = | + 1, (61) 
and A 2 as a function of the self-similar variable x/t is determined by the equation 

v 2 {v/2-l,X 2 ,X d _,v/2 + l) = - t . (62) 
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FIG. 7: Qualitative behavior of the Riemann invariants in the partial (attached) downstream dispersive shock. 



In the linear limit A2 — * Ai (i.e. m — > 0) we have 

v 2 (X 1 ,X 1 , A 3 , A 4 ) = A! + |(A 3 + A 4 ) + 2( V Al ! (A4 T Al) • ( 63 ) 

zAi — A3 — A4 

This velocity determines the speed si of the trailing edge of the downstream dispersive shock and it is not difficult 
to check using formulae (|43p . (|46[) that for V m -C 1 it is negative, that is at least for weak potentials the trailing edge 
cannot be located in the downstream region. Hence, in the equation (|62[) the variable A2 is limited from below by 
some cut-off value A|, and, as a result, the downstream shock gets attached to the barrier at its boundary x « I . 
Then one can find X 2 in the approximation described above by taking l/t — ► in Eq. (162j) , which yields 

v 2 {v/2-l,X 2l X d _,v/2 + l) = 0. (64) 

Since v 2 has the meaning of nonlinear group velocity, equation (|64[) can be interpreted as an asymptotic condition 
that the modulated wave "stops" at x — and this is why it can be directly connected (on the level of the Riemann 
invariants, i.e. in the "averaged" sense) to the stationary hydraulic transition, which on the modulation length scale 
can be viewed as a discontinuity located at x = 0. In more precise terms, relation (|64|) means that for the solution 
under study the characteristic of the Whitham equations dx/dt = v 2 coincides with the characteristic dx/dt = of the 
dispersionless equations (|3Tj) in the hydraulic approximation at x = so the modulation solution can be "terminated" 
at x = and the free-boundary matching conditions ()57|) at the (nonexistent) trailing edge x~ < can be replaced 
with the boundary conditions at x = 0: 

X 3 = Xi, A 4 =|u + 1, A x = iu-1 at x = 0. (65) 

The modulation solution is then considered in the upper right quarter of x, i-plane, x > 0, t > 0. One should mention 
that the similar situation occurs in the problem of the transcritical shallow water fluid flow past an obstacle described 
by the forced KdV equation where the "partial undular bore" attached to the obstacle is generated in the upstream 
flow [2^, [2^] (see also [13]). Qualitative behavior of the Riemann invariants in the attached downstream dispersive 
shock is shown in Fig. 7. Equation ([64]) determines A|, and therefore, via ([51]) . the value of the modulus to* at which 
the downstream dispersive shock is generated at x — 0. This value is equal to 



(A* -1„ + !)(!„ + l-Al) 

a v + i-x* 2 )(x d _-±v + iy 



(66) 



Within the partial dispersive shock the modulus changes in the interval m* < m < 1. The dependence of the cut-off 
modulus to* on the BEC flow velocity v is shown in Fig. 8a for V m = 0.5. One can see that for this case the cut-off 
modulus ranges in the interval < m* < 0.75 and for the transcritical flow velocities v sufficiently close to the lower 
boundary u_ one can treat the downstream dispersive shock as a dark soliton train slowly propagating to the right 
relative to the barrier. This is the dark soliton train that should remain within the finite elongated BEC for some time 
after the potential has been swept and, therefore, could be observed in the experiment. Therefore it is instructive to 
study its parameters in more detail. 

First we calculate the downstream soliton emission rate. This can be done using the formula 



/* = U*/L* , 



(67) 
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FIG. 8: a) Cut-off modulus m* in the downstream dispersive shock (dark soliton train) and b) Downstream dark soliton 
emission rate /* at the potential location vs oncoming BEC velocity v. Both figures correspond to V m = 0.5. At v ~ 1.4 the 
downstream dispersive shock wave detaches from the potential. Crosses in (b) correspond to the numerical simulations data. 
The dashed line in (b) corresponds to the frequency at the trailing edge of the detached downstream dispersive shock. 



where 

U* = |(Ai + A; + A 3 + A 4 ) = |0 + \t + A*) , 

L* = L(X x ,^XsM = ( 68 ) 

^(|i> + l-A$)(A, + 

are the phase velocity ([51]) and the wavelength (f54|) respectively, calculated at the soliton train generation point x = 0. 

Dependence f*(v) for V m — 0.5 is shown in Fig. 8b where it is also compared with the data from the direct 
numerical simulation using the potential (fl~3"]) with V m = 0.5, a = 2. Some features seen on this plot deserve 
additional explanation. Firstly we notice that the analytically obtained value for the soliton emission frequency is 
finite for the flow velocities close to (but slightly greater than) the lower transcritical boundary V- « 0.2. At the same 
time, the numerically obtained values of /* suggest that the rate of soliton production tends to zero as one approaches 
the lower transcritical region boundary i>_ ps 0.2 (we did not observe any soliton generation at v = 0.25 for the time 
range up to t — 50). The disagreement between the analytical solution and actual behaviour of the soliton emission 
frequency is due to the failure, close to the transcritical region boundary v — i>_, of our main assumption about the 
existence of the local transcritical hydraulic solution forming discontinuities with the equilibrium basic flow (see the 
discussion in Section Bl). We note that the detailed analysis of the dynamical scaling law for the soliton emission 
frequency for near-critical NLS flows through short-range (delta- function) potentials was performed in [3(|, where 
the frequency was shown to vanish as S 1 ^ 2 , S being the deviation of the controlling parameter (for instance, potential 
strength) from its critical value. 

The second "non-standard" feature in Fig. 8b is an abrupt change, back to zero, of the emission frequency /* for 
v w 1.4. This change does not constitute the cease of the soliton generation downstream but simply reflects the 
fact that the downstream dispersive shock gets detached from the obstacle potential for velocities greater than some 
v = v* , so there are no waves generated at the potential location at x = (where the frequency /* is defined). Indeed, 
as v increases within the transcritical region [t>_,t> + ], the speed si of the trailing edge (computed formally by (|63p) 
can change the sign from minus to plus (see Fig. 9), which implies that for some v = v* the downstream dispersive 
shock must detach from the barrier. The detachment threshold velocity is determined from the condition 

s d _(v*) = 0, (69) 

where s_(v) — V2{m = 0) is given by equation (f63|) . For V m — 0.5 this velocity is v* ps 1.4 (which can also be clearly 
seen in Figs. 8a, b). 

One of the physical consequences of the downstream dispersive shock detachment from the obstacle is the zero 
frequency of oscillations for the drag force (see Section 2D below). At the same time, one should note that at the 
point of the detachment, the amplitude of the dispersive shock vanishes, so the described discontinuity in the emission 
frequency at the obstacle has no practical significance. We stress that the solitons keep get generated downstream for 
v* < v < v+i but the generation point — the trailing edge of the downstream dispersive shock — now moves away from 
the potential with constant velocity si . The corresponding wave frequency at the moving generation point is shown 
in Fig. 8b by the dashed line. 
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FIG. 9: Dependence of the downstream dispersive shock trailing edge speed si vs transcritical BEC velocity v calculated by 
formula (|63[) . The detachment point v = v* « 1.4 is found from the condition si = 0. 
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FIG. 10: Leading soliton parameters: a) amplitude A d and b) speed s+ — in the downstream dispersive shock (dark soliton 
train) vs oncoming BEC velocity v for V m = 0.5. Solid line: modulation solution; Crosses: direct numerical solution. 



The limit A2 — > A3 (m — > 1) corresponds to the leading edge x + of the downstream dispersive shock. The amplitude 
of the leading soliton is given by 



and it moves with the velocity 



A d = (A 4 - A 3 )(A 3 -\ 1 ) = (2 + v- u d )(u d - v), 



s d + =v 2 (\ 00 1 \ d _,\ d _,\™) =u d -l. 



(70) 



(71) 



The dependencies A d (v) and si_(v) for V m = 0.5 are shown in Fig. 10 along with the corresponding numerical 
simulations data. In Figs. 11, 12 we plot the same set of quantities: A d , s 1 ^, m* and /* — as functions of the potential 
strength V m for fixed v — 1. One can see that our solution predicts the main physical parameters of the downstream 




(b) 

0.20 



1.0 \/ 



FIG. 11: a) Cut-off modulus m* in the downstream dispersive shock and b) Downstream dark soliton emission rate /* vs 
potential strength V m . Both figures correspond to v = 1.0. Crosses in (b) correspond to the numerical simulations data. 
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FIG. 12: Leading soliton parameters: a) amplitude A d and b) speed s+ — in the downstream dispersive shock (dark soliton 
train) vs potential strength V m - Both figures correspond to v — 1.0. Solid line: modulation solution; Crosses: direct numerical 
solution. 



wave quite well provided the flow velocity is not too close to the critical values v± and the potential strength V m is 
not too large. 

Summarizing, the downstream dispersive shock occupies the region 

si < x < 4 ■ t, (72) 
where si = for u_ < v < v* and si = si for v* < v < v + . 

2. Upstream dispersive shock 

The calculations in this case are very similar to those for the downstream dispersive shock. Now the constant 
Riemann invariants are equal to 

A 1 = A1, A 3 = A! 3 = v/2- 1, X 4 = X+=v/2 + l, (73) 

and A2 as a function of x/t is determined by the equation 

v 2 {Xt,X 2 ,v/2-l,v/2 + l)=x/t. (74) 

The zero-amplitude trailing edge with A2 — > Ai propagates with the velocity 

si = v 2 (X u _,X1,v/2-l,v/2 + l) = 2A1 - - 1 (75) 

A_ — v/2 

which is always negative. The limit \ 2 — ► A3 corresponds to the leading (soliton) edge. The amplitude of the leading 
soliton is equal to 

A u = 2(v - u u ), (76) 

and it moves with the velocity 

s l + = v 2 {X u _ , v/2 - 1, v/2 - 1, v/2 + 1) = i(u u + v-2). (77) 
In the case of of small V m C 1 we get on using expansion (|43p that 

s+ = K*-l)-/f, (78) 

which implies that in the transcritical interval (|44]) s+ positive. Thus there exists the range of velocities v for which 
the upstream shock is also attached to the barrier and realized only partially. In contrast to the attached downstream 
dispersive shock, in the partial upstream wave the modulus m varies between and some cut-off value m** < 1, i.e. 
this wave can be viewed as a nonlinear oscillatory tail rather than solitary wave train. 
The value A2 = A2*, is found from the equation 

« 3 (A*,^*,t;/2-l,«/2 + l)=0 1 (79) 
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and determines the value of the cut-off modulus m** via (|5ip. The threshold velocity v — v** at which the upstream 
dispersive shock gets attached to the potential is found from the condition s" = and is determined implicitly by 
the equation 

u u (v) = 2-v, (80) 

so that in the interval of velocities 

v** < v < v+ (81) 
the upstream dispersive shock is attached to the barrier. Thus, the upstream shock occupies the region 

si ■ t < x < s** • t, (82) 
where s+* = for the velocities v in the interval (f8Tj) and s** = s" for u_ < v < v**. 



D. Consolidated wave pattern 



Putting together the analytical results of this Section, we obtain full asymptotic description of the wave pattern 
generated in the transcritical ID flow of a BEC past a wide penetrable barrier. The pattern consists of two dispersive 
shocks propagating upstream and downstream the barrier and connected with each other by the transcritical hydraulic 
transition localized over the potential barrier spatial range. The behavior of the density and velocity within the 
dispersive shocks is obtained by the substitution of the slowly varying similarity modulation solutions (|61|) , (|62p and 
([73| , (|74|) for Xj into the rapidly oscillating traveling wave solution (|50|) and (l52|) . One should note that the matching 
conditions (|57p used in the construction of the modulation solution only guarantee the continuous matching of the 
average flow at the dispersive shock boundaries. To get exact matching of the rapidly oscillating wave field within 
the dispersive shocks with the constant (or slowly varying) flow outside, one should use the higher order analysis. 
The "weak limit" formulation of the dispersive shock problem used above is now well established owing to the studies 
based on the Lax-Levermore-Venakides rigorous approach (see [U, [42[ and references therein). 

The obtained combined modulated/hydraulic solution is shown in Fig. 13 at t = 30 for the potential with V m = 0.5 
and the oncoming flow velocity v = 1 (phase adjustments within the accuracy of the modulation theory are made to 
ensure continuity of the graph at the boundaries of the dispersive shocks). 

We have also performed direct numerical integration the Gross-Pitaevskii equation ([3]) with boundary conditions 
The simulations were performed using two different methods: the classical finite difference explicit scheme [iij 
and the quasi-spectral split step method [45j]. Both methods gave the same results. In Fig. 14 the numerical solution 
of the GP equation is plotted for the potential (|13[) with V m = 0.5 and a — 2. One can see excellent agreement 
between the wave patterns in Figs. 14 and 15. 

An additional small wave located about x — 60 in the numerical solution (Fig. 14) is due to the generation of a 
small-amplitude right-propagating wave packet formed from the Bogoliubov linear waves created by a switching on 
the obstacle's potential (see (43[). The group velocity v g — duj/dk of the Bogoliubov waves obeying the dispersion 
relation uj{k) — k\J\ + k 2 /4 is always greater than the sound speed c s = ^fp = 1 and tends to this value in the long 
wavelength limit k — + 0. Taking into account that this wave packet is convected by the flow with velocity w, we find 
that the wave packet occupies the region x > (c s + v) ■ t and its predicted position is about x > 60 for the chosen 
parameters with c s = v = 1 and t = 30 in agreement with the numerical results. Besides this wave packet, one 
can notice a tiny rarefaction wave right before the wave packet. It is formed due to the small discrepancy between 
the upstream and downstream values of the Riemann invariant A + (see Fig. 5 and the explanation of the closure 
conditions in subsection II. C). It is not difficult to show using standard hydrodynamic reasoning that to leading order 
the density jump across this rarefaction wave should be Ap ss S = (V m /6) 3 ^ 2 (see ([40|) ) while the rarefaction wave 
speed is calculated as u + ^fp ~ v + 1. For the parameters V m — 0.5, v — 1 used in our numerical simulation this 
implies Ap ps 0.025 and at t — 30 the predicted position of the rarefaction wave is about x = 60. Both predictions 
completely agree with the numerical solution. 

The agreement between the analytical and numerical solutions seen in Figs. 13 and 14 is especially remarkable in 
view of the relatively small width of the potential, a = 2, used in the numerical simulations. This width is comparable 
with the dispersion length in the system, which is of order of unity, so the formal requirement I — cr/2 3> 1 of the 
applicability of the local hydraulic solution is clearly violated. The robustness of the hydraulic solution here looks 
quite surprising and deserves special attention. In this regard we note that our analytical construction implies that 
two potentially conflicting requirements should be satisfied: the potential barrier should be (i) broad enough for the 
hydraulic approximation to be applicable but (ii) not too broad for the similarity modulation solution could be used 
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FIG. 13: Plot of the combined analytical (hydraulic + modulated oscillatory) solution for the condensate density distribution 
in the wave pattern generated by a BEC flow with v — 1 through the potential barrier with V m = 0.5 at t = 30. 



P 

1.6 
1.4 
1.2 
1 

0.8 
0.6 
0.4 
0.2 


-100 -80 -60 -40 -20 20 40 60 80 100 



FIG. 14: Numerical simulation of the condensate density distribution as a function of the space coordinate x in the wave 
pattern generated by a BEC flow with v = 1 through the potential barrier (113[) with Vm = 0.5, a = 2. The evolution time is 
equal to t = 30. 



for the description of the dispersive shocks (i.e. the characteristic time of the establishing of the steady transcritical 
hydraulic solution should be much less than the characteristic time of the formation of the dispersive shock) . While 
it might look that these requirements are difficult to satisfy simultaneously, our numerical simulations show that the 
resulting analytical solution works quite well when a = O(l). We have also performed numerical simulations for v = 1 
with potentials of different shapes with the conclusion that for potentials with V m < 1 and u = 0(1) the parameters 
of the oscillatory structure almost do not depend on the actual potential width and shape. This implies that for a 
reasonably broad range of the barrier potentials the whole wave pattern is characterized only by the potential strength 
V m and the flow velocity v, which agrees with the parametrization in our analytical solution. On the other hand, 
our simulations with very broad potentials a = O(10) show that the the steady hydraulic transition with constant 
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FIG. 15: Dependence of the velocities V- and v + at the boundaries of the transcritical region (solid lines) and of the velocities 
of 11* and «** (dashed lines) on the maximum value V m of potential. 



jumps outside the potential does not form within a finite time interval so the developed quantitative description of 
the dispersive shocks with the aid the similarity modulation solutions does not apply. 

As we see in Figs. 13 and 14, the parameters v = 1, V m — 0.5 correspond to the case when the downstream 
dispersive shock is attached to the obstacle. It is natural to ask at which values of the parameters v, V m the upstream 
dispersive shock gets detached and whether there exists the region of the parameters when both shocks are detached 
from the obstacle. To answer these questions, we first notice that the downstream dispersive shock detaches from 
the obstacle when the velocity of the trailing (m = 0) edge of the shock v 2 {v/2 — l,v/2— 1, A_, v/2 + 1) defined by 
Eq. (163|) vanishes. This condition gives the equation 

\v 2 - (9 + \ d _)v - (X d _) 2 + 6Xt + 11 = 0. (83) 

Taking into account Eq. (l4Tj) we find Al = u d /2 — = u d — v/2 — 1, and cast this equation to the form 

v 2 - 12v - (u d ) 2 + 8u d + 4 = 0, (84) 

which gives the value v* of the "detachment" velocity 

v* = 6 - ^(^-4)2 + 16, (85) 
where u d is the greater root of the equation (I42[) . In the weak potential limit V m C 1 we obtain the series expansion 

V =1+ 2\j — - — + ---- (86) 

The upstream shock detaches from the obstacle at velocity v** which satisfies the equation (|50|) and again in the weak 
potential limit we obtain 

*"-i+5i/?-ir + -- < 87 > 

It is easy to see that the region v* < v < v** is located inside the transcritical region (14"4"1) and is relatively narrow; 
its position is illustrated in Fig. 15. 

We have verified the above predictions by constructing numerical solutions of the GP equation with V m = 0.5 and 
v = 1.4 (v* < v < v**), v — 1.8 (v > v**). The respective plots are presented in Fig. 16. 

In Fig. 16a both dispersive shocks are completely developed and actually they both are detached from the obstacle 
although because of the small difference between the soliton edge velocity of the upstream shock and the velocity of the 
small amplitude edge of the downstream shock, it takes very long time to reach a well developed hydraulic transition 
solution located near x — solution. Nevertheless we see that both dispersive shocks are not cut off at the edges close 
to the obstacle which means their detachment from the obstacle. In contrast, in Fig. 16b the downstream dispersive 
shock is detached from the obstacle whereas the upstream shock is attached and is cut off towards the soliton edge. 
However, one should notice that the parameters of the dispersive shocks for the velocity v near the upper boundary 
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v + of the transcritical region do not agree well enough with the analytical predictions. This disagreement has already 
been discussed in Section Bl and is due to the violation of our main supposition that the flow forms a steady hydraulic 
transition over the potential range interval with the jumps at the both its sides: here we cannot neglect the time of 
forming of the hydraulic solution compared with the time of the development of the dispersive shocks and therefore, 
the self-similar solutions used in the analytical theory are not accurate enough. In spite of this reservation, the 
developed theory qualitatively agrees with numerics even in this region: one can easily distinguish in the numerically 
obtained pattern all the characteristic ingredients of our analytical construction, namely, the smooth transition region 
over the potential range and downstream and the upstream dispersive shocks. 



E. Drag force 

We now consider the drag force, i.e. the force exerted on the BEC due to its motion through the potential barrier. 
This force can be calculated as the spatial mean value of the operator dV(x)/dx over the condensate wave function 
(see [1^ for a detailed derivation), 

oo oc 

f dV f dV 

Fdrag = / tp —t()dx = / p — dx. (88) 

— oo — oo 

For subcritical and supercritical flows, p = p(x) is given by the hydraulic solution (jT^J) which connects smoothly to 
p = 1 as |x| — > oo and satisfies the system ([5]). Then integrating the associated Bernoulli equation 

(pu 2 + p 2 /2) x + P V X =0 (89) 

from — oo to +00, and using that p — > 1, u — * v as \x\ — > 00 we obtain Fdrag = which is the expression of BEC 
superfluidity at sub- and supercritical velocities. Strictly speaking, the superfluidity is exact for v < V- < c s = 1 
where no excitations are generated. In the supersonic region v > v+ > c s = 1 the generation of excitations exists but 
it is exponentially small for a slowly varying obstacles 's potential [251 ] and this generation is neglected in the hydraulic 
approximation used here. 

For the transcritical regime, V- < v < v+, there is no global hydraulic solution and the integral (1551) can be split 
into three parts 

— l l CO 

f dV f dV f dV 

F d rag= / Pd{x,t)—dx+ p tr (x) — dx+ p u (x,t) — dx. (90) 

-oo -i I 

Here p u ,d(x,t) are the unsteady upstream and downstream solutions describing the density behavior in respective 
dispersive shocks and ptr{x) is the local transcritical hydraulic solution (|30|) defined on the interval (— I, I). Assuming 
that the potential V(x) sufficiently rapidly decays for |x| > I together with its first derivative V x , one can neglect the 
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FIG. 17: The drag force Fdrag vs: a ) BEC oncoming flow v in the transcritical region (v-,v+) and b) potential strength V m 
for fixed v = 1. Solid lines are drawn according to the approximate Eq. (|97[) and dots correspond to the full numerical solution 
of Eq. (PI). 



contributions of the first and third integrals in (|90[) . The remaining second integral can be evaluated again with the 
use of the Bernoulli equation ([89]) to obtain 

F drag = p"K) 2 + \{P U ? - P d {u d f ~ |(/) 2 . (91) 

If the downstream or upstream shock is attached to the hydraulic solution, then the limiting values of the density 
and the flow velocity in the right-hand side of Eq. (f9"Tj) oscillate with time according to the dispersive shock solution 
(|50|> . (|53j) considered at x — and the parameter m equal to m* or m** depending on whether the downstream or 
upstream dispersive shock is attached. The frequency of the drag force oscillations for the downstream attachment 
case is given by formula (|67p (a similar formula can be easily obtained for the upstream attachment case as well). 
Averaging of the expression (|9ip over time yields the mean value of the drag force. The situation simplifies greatly 
when both dispersive shocks are detached from the obstacle and p u,d and u u ' d are given by the limiting values of the 
hydraulic solution. Although the corresponding region of the potential maximum V m and velocity v values is rather 
narrow, the discussion of this case is quite instructive and enables one to estimate the accuracy of the drag force series 

1/2 

expansion in powers of V m for small V m ■ 

As was shown above, w = u u,d are two roots of the equation (|42p which can be rewritten in a more convenient form 
after introduction of new variables 

w = u {l-z), P = ul{l + z/2f 1 u Q = {v + 2)/3, e = V m /v% = 9V m /(v + 2) 2 , (92) 

so that Eq. (|42p takes the form 

|(1 - zf + (1 + z/2) 2 - |(1 - z) 2 / 3 (l + z/2) 4 / 3 = e. (93) 
One can easily derive series expansions of the roots of this equation in powers of e 1 / 2 to obtain 

h . . . , z d = -\ V ... , 94) 

3 18 ' V 3 18 ' V ' 

where the first order terms reproduce actually Eqs. I|43p after taking into account inequalities (144p . With the use of 
Eq. Q we represent Eq. (gTJ) as 

Fdrag = 4 [G(z U ) - G{z d )\ (95) 

where 

G(z) = (1 - z) 2 (l + z/2) 2 + |(1 + z/2) 4 . (96) 
Then substitution of (|94p into ([95]) yields with the accepted accuracy the expression 

(97) 




= £ 5/2 



6V6 / V 6 



v \ 3 / 2 



A(v + 2) - \v n 
6 



where v varies in the transcritical region (144p . In Fig. 17 we compare the plots of dependence of Fdrag on v and V m 
according to Eq. (f9"T|) and calculated by means of the exact numerical solution of Eq. (l9"3")) . As we see, the accuracy 
is very good even for V m — 1. 
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III. DISCUSSION 

In the experiment [2lT | it was found that the solitons are generated by a moving potential barrier in the interval of 
velocities 

0.3mm/s < v < 0.9mm/s. (98) 

This result agrees qualitatively with existence of the finite interval (|14[) for which the expanding dispersive shocks are 
generated. However, we encounter a quantitative contradiction if we accept the value of the sound velocity calculated 
in [2ll ] c s = 2.1mm/s as correct, because in our non-dimensional units the sound velocity is equal to unity and 
hence it must be located inside the interval (fT4")l . (see Fig. 2), or, in dimensional units, inside the interval (pS)) . This 
disagreement can be explained by noticing that the above value of the sound velocity was calculated in [2l| according 
to the expression 

(99) 

where po is the condensate density at the center of the trap. Here 

g = 4irh 2 a s /m (100) 

is the effective coupling constant in the BEC consisting of atoms with mass m and s-wave scattering length a s . But 
this expression is correct only for a rarefied enough condensate confined to the cigar-shaped trap with "frozen" radial 
motion and this condition was not fulfilled in [21| . 

Dynamics of a dense BEC is described by the full 3D GP equation, which can be reduced to some effectively ID 
systems much more complicated than the NLS equation (fT]). For example, the variational approach to the dynamics 
of the dense BEC was developed in [46| where it was shown that the sound velocity along the axial direction of the 
trap is given by the expression (see Eqs. (67) and (70) in [46| ) 

Cl = c ».(l±Mgf, {101) 



(1 + 2Gf/ 4 



where the parameter G is calculated by the formula 




to 








Uv 





G = tk\\ 1+ \tk) +tk\- ( 102 ) 



Here a± = ^h/moj± is the radial "oscillator length" and £ = hi \/2mpy,q is the healing length. Eq. (|101[) reduces to 
Eq. (|99|) in the limit G <C 1 of a rarefied BEC. In the experiment [2J| these parameters were equal to a± = 0.73- 10~ 4 cm 
and £ = 0.17 • 10~ 4 cm, so that Eq. (|102[) gives G — 11.4, that is the experiment [2l| corresponds to the opposite 
limit of dense BEC. The sound velocity calculated according to Eq. (|101|) is equal to c s — 0.8mm/s, and this value 
agrees much better with the interval (|98p . Thus, for quantitative description of the experiment [2l| the theory of 
dispersive shocks in a dense BEC should be developed what is beyond the scope of this paper. Therefore we shall 
confine ourselves here to some qualitative remarks only. 

The famous Landau criterion Q for the loss of superfluidity was based on the consideration of linear excitations 
only and it contradicted to the experiments with liquid Hell. This discrepancy was explained by Feynman by taking 
into the consideration the formation of such nonlinear structures as vortices in the flow of a superfluid in capillaries or 
past obstacles. However, the notion of the threshold velocity below which the flow is superfluid has not been changed 
by this modification of the theory. Taking into account the generation of dispersive shocks in 2D situation 0, [l(| 
did not change this notion either, since the stationary spatial dispersive shocks are generated by the supersonic flow 
of BEC only. The results of the present paper, as well as of the previous papers [23|, [H, [2^, [2(| |27[, show that the 
situation can be more subtle in the case of ID flow. In this case, the flow past a wide barrier leads to the generation of 
dispersive shocks in a finite interval of the flow velocities bounded not only from below but also from above. Moreover, 
the lower boundary value of the velocity V- can become equal to zero for strong enough barriers, that is even very slow 
motions could lead to the generation of solitons. This observation is in striking contrast with the standard reasonings 
based on the linear theory of excitations. 
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